'data.frame': 544 obs. of 4 variables:
$ height: num 152 140 137 157 145 ...
$ weight: num 47.8 36.5 31.9 53 41.3 ...
$ age : num 63 63 65 41 51 35 32 27 19 54 ...
$ male : int 1 0 0 1 0 1 0 1 0 1 ...
Un cours en R et Stan avec brms
Cours n°01 : Introduction à l’inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux (généralisés)
Cours n°09 : Examen final
\[\newcommand\given[1][]{\:#1\vert\:}\]
\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i}&= \alpha + \beta x_{i} \\ \alpha &\sim \mathrm{Normal}(60, 10) \\ \beta &\sim \mathrm{Normal}(0, 10) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 1) \end{align} \]
Objectif de la séance : comprendre ce type de modèle.
Les constituants de nos modèles seront toujours les mêmes et nous suivrons les trois mêmes étapes :
'data.frame': 544 obs. of 4 variables:
$ height: num 152 140 137 157 145 ...
$ weight: num 47.8 36.5 31.9 53 41.3 ...
$ age : num 63 63 65 41 51 35 32 27 19 54 ...
$ male : int 1 0 0 1 0 1 0 1 0 1 ...
\[h_{i} \sim \mathrm{Normal}(\mu, \sigma)\]
\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]
Contraintes : Certaines valeurs soient fortement probables (autour de la moyenne \(\mu\)). Plus on s’éloigne, moins les valeurs sont probables (en suivant une décroissance exponentielle).
\[ y = \exp \big[-x^{2} \big] \]
On étend notre fonction aux valeurs négatives.
\[ y = \exp \big[-x^{2} \big] \]
Les points d’inflection nous donnent une bonne indication de là où la plupart des valeurs se trouvent (i.e., entre les points d’inflection). Les pics de la dérivée nous montrent les points d’inflection.
\[ y = \exp \bigg [- \frac{1}{2} x^{2} \bigg] \]
Ensuite on standardise la distribution de manière à ce que les deux points d’inflection se trouvent à \(x = -1\) et \(x = 1\).
\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} x^{2} \bigg] \]
On insère un paramètre \(\sigma^{2}\) pour contrôler la distance entre les points d’inflection.
\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]
On insère ensuite un paramètre \(\mu\) afin de pouvoir contrôler la position (la tendance centrale) de la distribution.
\[ y = \frac{1}{\sqrt{2 \pi \color{steelblue}{\sigma^{2}}}} \exp \bigg[-\frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]
Mais… cette distribution n’intègre pas à 1. On divise donc par une constante de normalisation (la partie gauche), afin d’obtenir une distribution de probabilité.
Nous allons construire un modèle de régression, mais avant d’ajouter un prédicteur, essayons de modéliser la distribution des tailles.
On cherche à savoir quel est le modèle (la distribution) qui décrit le mieux la répartition des tailles. On va donc explorer toutes les combinaisons possibles de \(\mu\) et \(\sigma\) et les classer par leurs probabilités respectives.
Notre but, une fois encore, est de décrire la distribution postérieure, qui sera donc d’une certaine manière une distribution de distributions.
On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).
\[\color{steelblue}{\mu \sim \mathrm{Normal}(178, 20)}\]
On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).
\[\color{steelblue}{\sigma \sim \mathrm{Uniform}(0, 50)}\]
library(ks)
sample_mu <- rnorm(1e4, 178, 20) # prior on mu
sample_sigma <- runif(1e4, 0, 50) # prior on sigma
prior <- data.frame(cbind(sample_mu, sample_sigma) ) # multivariate prior
H.scv <- Hscv(x = prior, verbose = TRUE)
fhat_prior <- kde(x = prior, H = H.scv, compute.cont = TRUE)
plot(
fhat_prior, display = "persp", col = "steelblue", border = NA,
xlab = "\nmu", ylab = "\nsigma", zlab = "\n\np(mu, sigma)",
shade = 0.8, phi = 30, ticktype = "detailed",
cex.lab = 1.2, family = "Helvetica")[1] 162.8648
On veut calculer la probabilité d’observer une certaine valeur de taille, sachant certaines valeurs de \(\mu\) et \(\sigma\), c’est à dire :
\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]
\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]
Ou avec une fonction maison…
\[ \color{purple}{p(\mu, \sigma \given h)} = \frac{\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i} \given \mu, \sigma)\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50) \mathrm{d} \mu \mathrm{d} \sigma}} \]
\[ \color{purple}{p(\mu, \sigma \given h)} \propto \prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)} \]
Il s’agit de la même formule vue lors des cours 1 et 2, mais cette fois en considérant qu’il existe plusieurs observations de taille (\(h_{i}\)), et deux paramètres à estimer : \(\mu\) et \(\sigma\).
Pour calculer la vraisemblance marginale (en vert), il faut donc intégrer sur deux paramètres : \(\mu\) et \(\sigma\). On réalise ici encore que la probabilité a posteriori est proportionnelle au produit de la vraisemblance et du prior.
# on définit une grille de valeurs possibles pour mu et sigma
mu.list <- seq(from = 140, to = 160, length.out = 200)
sigma.list <- seq(from = 4, to = 9, length.out = 200)
# on étend la grille à toutes les combinaisons possibles de mu et sigma
post <- expand.grid(mu = mu.list, sigma = sigma.list)
# calcul de la log-vraisemblance (pour chaque combinaison de mu et sigma)
post$LL <-
sapply(
1:nrow(post),
function(i) sum(dnorm(
d2$height,
mean = post$mu[i],
sd = post$sigma[i],
log = TRUE) )
)
# calcul de la probabilité a posteriori (non normalisée)
post$prod <-
post$LL +
dnorm(x = post$mu, mean = 178, sd = 20, log = TRUE) +
dunif(x = post$sigma, min = 0, max = 50, log = TRUE)
# on "annule" le log et on standardise par la valeur maximale (pour éviter les erreurs d'arrondi)
post$prob <- exp(post$prod - max(post$prod) ) mu sigma LL prod prob
1 140.8040 4.577889 -2958.659 -2968.215 0.000000e+00
2 140.1005 8.070352 -1787.919 -1797.541 4.138624e-248
3 149.8492 8.472362 -1277.453 -1286.270 4.560547e-26
4 140.9045 6.286432 -2071.760 -2081.307 0.000000e+00
5 158.3920 8.246231 -1258.081 -1266.388 1.967155e-17
6 140.0000 8.195980 -1778.828 -1788.460 3.638581e-244
7 143.8191 7.366834 -1596.980 -1606.267 4.852820e-165
8 146.4322 4.678392 -1883.309 -1892.381 2.682224e-289
9 154.2714 5.231156 -1291.017 -1299.547 7.815693e-32
10 149.0452 6.663317 -1350.199 -1359.074 1.098510e-57
11 150.2513 5.557789 -1375.406 -1384.195 1.351857e-68
12 149.6482 7.819095 -1289.953 -1298.784 1.675922e-31
13 144.4221 7.291457 -1563.391 -1572.627 1.976531e-150
14 155.9799 6.989950 -1230.123 -1238.556 2.404888e-05
15 156.8844 7.266332 -1238.257 -1246.641 7.409259e-09
16 142.0101 7.668342 -1693.624 -1703.070 4.418555e-207
17 152.0603 4.653266 -1402.852 -1411.519 1.836562e-80
18 140.4020 6.261307 -2142.119 -2151.712 0.000000e+00
19 154.7739 6.939698 -1223.939 -1232.440 1.089168e-02
20 149.1457 8.246231 -1297.723 -1306.590 6.826903e-35
Under the hood : Stan est un langage de programmation probabiliste écrit en C++, et qui implémente plusieurs algorithmes de MCMC : HMC, NUTS, L-BFGS…
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu;
real<lower=0> tau;
real eta[J];
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}Le package brms (Bürkner, 2017) permet de fitter des modèles multi-niveaux (ou pas) linéaires (ou pas) bayésiens en Stan mais en utilisant la syntaxe de lme4.
Par exemple, le modèle suivant :
\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \alpha_{\text{subject}[i]} + \alpha_{\text{item}[i]} + \beta x_{i} \\ \end{align} \]
Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.
La partie gauche représente notre variable dépendante (ou outcome, i.e., ce qu’on essaye de prédire). Le package brms permet également de fitter des modèles multivariés (plusieurs outcomes) en les combinant avec mvbind().
Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.
Par défaut brms postule une vraisemblance gaussienne. Ce postulat peut être changé facilement en spécifiant la vraisemblance souhaitée via l’argument family.
Lisez la documentation (c’est très enthousiasmant à lire) accessible via ?brm.
# générer le code du modèle en Stan
make_stancode(formula, ...)
stancode(fit)
# définir les priors
get_prior(formula, ...)
set_prior(prior, ...)
# récupérer les prédictions du modèle
fitted(fit, ...)
predict(fit, ...)
conditional_effects(fit, ...)
# posterior predictive checking
pp_check(fit, ...)
# comparaison de modèles
loo(fit1, fit2, ...)
bayes_factor(fit1, fit2, ...)
model_weights(fit1, fit2, ...)
# test d'hypothèse
hypothesis(fit, hypothesis, ...) Estimate Est.Error Q2.5 Q97.5
b_Intercept 154.6044 0.4212561 153.790605 155.416941
sigma 7.7610 0.3050007 7.204519 8.370314
Ces données représentent les distributions marginales de chaque paramètre. En d’autres termes, la probabilité de chaque valeur de \(\mu\), moyennée sur toutes les valeurs possible de \(\sigma\), est décrite par une distribution gaussienne avec une moyenne de \(154.6\) et un écart type de \(0.42\). L’intervalle de crédibilité (\(\neq\) intervalle de confiance) nous indique les 95% valeurs de \(\mu\) ou \(\sigma\) les plus probables (sachant les données et les priors).
Par défaut brms utilise un prior très peu informatif centré sur la valeur moyenne de la variable mesurée. On peut donc affiner l’estimation réalisée par ce modèle en utilisant nos connaissances sur la distribution habituelle des tailles chez les humains.
La fonction get_prior() permet de visualiser une liste des priors par défaut ainsi que de tous les priors qu’on peut spécifier, sachant une certaine formule (i.e., une manière d’écrire notre modèle) et un jeu de données.
Family: gaussian
Links: mu = identity
Formula: height ~ 1
Data: d2 (Number of observations: 352)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 154.60 0.42 153.80 155.42 1.00 3285 2855
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 7.77 0.30 7.21 8.37 1.00 2890 2245
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Family: gaussian
Links: mu = identity
Formula: height ~ 1
Data: d2 (Number of observations: 352)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 177.86 0.10 177.67 178.07 1.00 2985 2622
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 24.60 0.92 22.90 26.49 1.00 3157 2848
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
On remarque que la valeur estimée pour \(\mu\) n’a presque pas “bougée” du prior…mais on remarque également que la valeur estimée pour \(\sigma\) a largement augmentée. Nous avons dit au modèle que nous étions assez certain de notre valeur de \(\mu\), le modèle s’est ensuite “adapté”, ce qui explique la valeur de \(\sigma\)…
Le prior peut généralement être considéré comme un posterior obtenu sur des données antérieures.
On sait que le \(\sigma\) d’un posterior gaussien nous est donné par la formule :
\[\sigma_{\text{post}} = 1\ /\ \sqrt{n}\]
Qui implique une quantité de données \(n = 1\ /\ \sigma^2_{\text{post}}\). Notre prior avait un \(\sigma = 0.1\), ce qui donne \(n = 1\ /\ 0.1^2 = 100\).
On peut donc considérer que le prior \(\mu \sim \mathrm{Normal}(178, 0.1)\) est équivalent au cas dans lequel nous aurions observé \(100\) tailles de moyenne \(178\).
# A draws_df: 6 iterations, 1 chains, and 6 variables
b_Intercept sigma Intercept lprior lp__ density
1 155 8.0 155 -9.3 -1227 0.851
2 154 8.0 154 -9.3 -1227 0.891
3 155 7.3 155 -9.3 -1228 0.391
4 154 7.3 154 -9.3 -1229 0.122
5 154 7.2 154 -9.3 -1229 0.128
6 154 7.2 154 -9.3 -1231 0.033
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
Comment est-ce que la taille co-varie avec le poids ?
\[ \begin{align} h_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \beta x_{i} \\ \end{align} \]
mean sd 2.5% 97.5%
(Intercept) 113.8793936 1.91106523 110.1337746 117.6250126
weight 0.9050291 0.04204752 0.8226175 0.9874407
On considère un modèle de régression linéaire avec un seul prédicteur, une pente, un intercept, et des résidus distribués selon une loi normale. La notation :
\[ h_{i} = \alpha + \beta x_{i} + \epsilon_{i} \quad \text{avec} \quad \epsilon_{i} \sim \mathrm{Normal}(0, \sigma) \]
est équivalente à :
\[ h_{i} - (\alpha + \beta x_{i}) \sim \mathrm{Normal}(0, \sigma) \]
et si on réduit encore un peu :
\[ h_{i} \sim \mathrm{Normal}(\alpha + \beta x_{i}, \sigma). \]
Les notations ci-dessus sont équivalentes, mais la dernière est plus flexible, et nous permettra par la suite de l’étendre plus simplement aux modèles multi-niveaux.
\[ \begin{aligned} \color{orangered}{h_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i},\sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta x_{i}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(178, 20)} \\ \color{steelblue}{\beta} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]
Dans ce modèle \(\mu\) n’est plus un paramètre à estimer (car \(\mu\) est déterminé par \(\alpha\) et \(\beta\)). À la place, nous allons estimer \(\alpha\) et \(\beta\).
Rappels : \(\alpha\) est l’intercept, c’est à dire la taille attendue, lorsque le poids est égal à \(0\). \(\beta\) est la pente, c’est à dire le changement de taille attendu quand le poids augmente d’une unité.
Estimate Est.Error Q2.5 Q97.5
b_Intercept 113.8765923 1.92009947 110.188131 117.6504948
b_weight 0.9049581 0.04206323 0.820382 0.9846767
sigma 5.1016861 0.19447071 4.736904 5.4935198
Intercept 154.5910945 0.27375026 154.054875 155.1235651
lprior -12.4815498 0.01618103 -12.513986 -12.4502345
lp__ -1083.3808160 1.26795384 -1086.812949 -1081.9551330
\(\alpha = 113.89, 95\% \ \text{CrI} \ [110.15, 117.58]\) représente la taille moyenne quand le poids est égal à 0kg…
\(\beta = 0.90, 95\% \ \text{CrI} \ [0.82, 0.99]\) nous indique qu’une augmentation de 1kg entraîne une augmentation de 0.90cm.
Estimate Est.Error Q2.5 Q97.5
Intercept 154.6015921 0.27311256 154.0899441 155.1354267
weight.c 0.9051041 0.04201935 0.8230953 0.9874308
Après avoir centré la réponse, l’intercept représente désormais la valeur attendue de taille (en cm) lorsque le poids est à sa valeur moyenne.
d2 %>%
ggplot(aes(x = weight, y = height) ) +
geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
# geom_abline(intercept = fixef(mod4)[1], slope = fixef(mod4)[2], lwd = 1)
geom_smooth(
data = data.frame(fitted(object = mod4) ) %>% bind_cols(data.frame(weight = d2$weight) ),
aes(y = Estimate), stat = "identity",
se = FALSE, color = "black"
)# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 25, to = 70, by = 1) )
# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(fitted(mod4, newdata = weight.seq) ) %>% bind_cols(weight.seq)
# on affiche les 10 premières lignes de mu
head(mu, 10) Estimate Est.Error Q2.5 Q97.5 weight
1 136.5005 0.8919314 134.7632 138.2530 25
2 137.4055 0.8519943 135.7486 139.0836 26
3 138.3105 0.8122718 136.7299 139.9143 27
4 139.2154 0.7727971 137.7124 140.7353 28
5 140.1204 0.7336101 138.6924 141.5628 29
6 141.0253 0.6947596 139.6694 142.3904 30
7 141.9303 0.6563053 140.6504 143.2091 31
8 142.8353 0.6183211 141.6363 144.0411 32
9 143.7402 0.5808992 142.6096 144.8678 33
10 144.6452 0.5441558 143.5843 145.7049 34
Pour rappel, voici notre modèle : \(h_{i} \sim \mathrm{Normal}(\alpha + \beta x_{i}, \sigma)\). Pour l’instant, on a seulement représenté les prédictions pour \(\mu\). Comment incorporer \(\sigma\) dans nos prédictions ?
# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 25, to = 70, by = 1) )
# on récupère les prédictions du modèle pour ces valeurs de poids
pred_height <- data.frame(predict(mod4, newdata = weight.seq) ) %>% bind_cols(weight.seq)
# on affiche les 10 premières lignes de pred_height
head(pred_height, 10) Estimate Est.Error Q2.5 Q97.5 weight
1 136.4612 5.086325 126.3412 146.2835 25
2 137.4636 5.167476 127.2898 147.5171 26
3 138.2198 5.198687 128.1563 148.2731 27
4 139.3200 5.174066 129.1870 149.5011 28
5 140.1321 5.243154 129.7430 150.2355 29
6 141.0754 5.094878 130.9877 150.9205 30
7 141.8453 5.130618 131.8189 151.8517 31
8 142.8526 5.130884 133.0214 152.9146 32
9 143.6994 5.234287 133.4353 154.0481 33
10 144.5977 5.225459 134.2982 155.1205 34
d2 %>%
ggplot(aes(x = weight, y = height) ) +
geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
geom_ribbon(
data = pred_height, aes(x = weight, ymin = Q2.5, ymax = Q97.5),
alpha = 0.2, inherit.aes = FALSE
) +
geom_smooth(
data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity", color = "black", alpha = 0.8, size = 1
)Le paquet brms propose aussi les fonctions posterior_epred(), posterior_linpred(), et posterior_predict(), qui permettent de générer des prédictions à partir de modèles fittés avc brms. Andrew Heiss décrit de manière détaillée le fonctionnement de ces fonction dans cet article de blog.
Deux sources d’incertitude dans le modèle : incertitude concernant l’estimation de la valeur des paramètres mais également concernant le processus d’échantillonnage.
Incertitude épistémique : La distribution a posteriori ordonne toutes les combinaisons possibles des valeurs des paramètres selon leurs plausibilités relatives.
Incertitude aléatoire : La distribution des données simulées est elle, une distribution qui contient de l’incertitude liée à un processus d’échantillonnage (i.e., générer des données à partir d’une gaussienne).
Voir aussi ce court article par O’Hagan (2004).
[1] 3.346996e-17 1.000000e+00
Pourquoi standardiser les prédicteurs ?
Interprétation. Permet de comparer les coefficients de plusieurs prédicteurs. Un changement d’un écart-type du prédicteur correspond à un changement d’un écart-type sur la réponse (si la réponse est aussi standardisée).
Fitting. Quand les prédicteurs contiennent de grandes valeurs (ou des valeurs trop différentes les unes des autres), cela peut poser des problèmes de convergence (cf. Cours n°05).
\[ \begin{aligned} \color{orangered}{h_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(156, 100)} \\ \color{steelblue}{\beta_{1}, \beta_{2}} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]
À vous de construire et fitter ce modèle en utilisant brms::brm().
Family: gaussian
Links: mu = identity
Formula: height ~ 1 + weight.s + I(weight.s^2)
Data: d (Number of observations: 544)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 146.66 0.38 145.91 147.38 1.00 3621 2764
weight.s 21.40 0.29 20.83 21.96 1.00 3181 2802
Iweight.sE2 -8.41 0.28 -8.97 -7.87 1.00 3590 3289
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 5.78 0.17 5.46 6.13 1.00 3876 2793
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight.s = seq(from = -2.5, to = 2.5, length.out = 50) )
# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(fitted(mod6, newdata = weight.seq) ) %>% bind_cols(weight.seq)
pred_height <- data.frame(predict(mod6, newdata = weight.seq) ) %>% bind_cols(weight.seq)
# on affiche les 10 premières lignes de pred_height
head(pred_height, 10) Estimate Est.Error Q2.5 Q97.5 weight.s
1 40.53360 5.994749 28.98851 52.16675 -2.500000
2 46.75169 5.800168 35.68154 58.15220 -2.397959
3 52.98804 5.914122 41.60599 65.08939 -2.295918
4 59.31354 5.772298 48.18958 70.66286 -2.193878
5 64.94578 5.911210 53.32341 76.66299 -2.091837
6 70.76235 5.747479 59.35056 81.63207 -1.989796
7 76.32943 5.865234 64.81440 87.91037 -1.887755
8 81.61468 5.813679 70.28378 92.81092 -1.785714
9 86.88819 5.646336 75.66007 97.81014 -1.683673
10 91.69154 5.847410 80.35947 103.24958 -1.581633
d %>%
ggplot(aes(x = weight.s, y = height) ) +
geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
geom_ribbon(
data = pred_height, aes(x = weight.s, ymin = Q2.5, ymax = Q97.5),
alpha = 0.2, inherit.aes = FALSE
) +
geom_smooth(
data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity", color = "black", alpha = 0.8, size = 1
)Plusieurs méthodes pour calculer les tailles d’effet dans les modèles bayésiens. Gelman & Pardoe (2006) proposent une méthode pour calculer un \(R^{2}\) basé sur l’échantillon.
Marsman & Wagenmakers (2017) et Marsman et al. (2019) généralisent des méthodes existantes pour calculer un \(\rho^{2}\) pour les designs de type ANOVA (i.e., avec prédicteurs catégoriels), qui représente une estimation de la taille d’effet dans la population (et non basée sur l’échantillon).
Similar to most of the ES measures that have been proposed for the ANOVA model, the squared multiple correlation coefficient \(\rho^{2}\) […] is a so-called proportional reduction in error measure (PRE). In general, a PRE measure expresses the proportion of the variance in an outcome \(y\) that is attributed to the independent variables \(x\) (Marsman et al., 2019).
\[ \begin{aligned} \rho^{2} &= \dfrac{\sum_{i = 1}^{n} \pi_{i}(\beta_{i} - \beta)^{2}}{\sigma^{2} + \sum_{i=1}^{n} \pi_{i}(\beta_{i} - \beta)^{2}} \\ \rho^{2} &= \dfrac{ \frac{1}{n} \sum_{i=1}^{n} \beta_{i}^{2}}{\sigma^{2} + \frac{1}{n} \sum_{i = 1}^{n} \beta_{i}^{2}} \\ \rho^{2} &= \dfrac{\beta^{2} \tau^{2}}{\sigma^{2} + \beta^{2} \tau^{2}}\\ \end{aligned} \]
On a présenté un nouveau modèle à deux puis trois paramètres : le modèle gaussien, puis la régression linéaire gaussienne, permettant de mettre en relation deux variables continues.
Comme précédemment, le théorème de Bayes est utilisé pour mettre à jour nos connaissances a priori quant à la valeur des paramètres en une connaissance a posteriori, synthèse entre nos priors et l’information contenue dans les données.
La package brms permet de fitter toutes sortes de modèles avec une syntaxe similaire à celle utilisée par lm().
La fonction fitted() permet de récupérer les prédictions d’un modèle fitté avec brms.
La fonction predict() permet de simuler des données à partir d’un modèle fitté avec brms.
Sélectionner toutes les lignes du jeu de données howell correspondant à des individus mineurs (age < 18). Cela devrait résulter en une dataframe de 192 lignes.
Fitter un modèle de régression linéaire en utilisant la fonction brms::brm(). Reporter et interpréter les estimations de ce modèle. Pour une augmentation de 10 unités de weight, quelle augmentation de taille (height) le modèle prédit-il ?
Faire un plot des données brutes avec le poids sur l’axe des abscisses et la taille sur l’axe des ordonnées. Surimposer la droite de régression du modèle et un intervalle de crédibilité à 89% pour la moyenne. Ajouter un intervalle de crédibilité à 89% pour les tailles prédites.
Que pensez-vous du “fit” du modèle ? Quelles conditions d’application du modèle seriez-vous prêt.e.s à changer, afin d’améliorer le modèle ?
Imaginons que vous ayez consulté une collègue experte en allométrie (i.e., les phénomènes de croissance différentielle d’organes) et que cette dernière vous explique que ça ne fait aucun sens de modéliser la relation entre le poids et la taille… alors qu’on sait que c’est le logarithme du poids qui est relié (linéairement) à la taille !
Modéliser alors la relation entre la taille (cm) et le log du poids (log-kg). Utiliser la dataframe howell en entier (les 544 lignes). Fitter le modèle suivant en utilisant brms::brm().
\[ \begin{align*} &\color{orangered}{h_{i} \sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ &\mu_{i}= \alpha + \beta \cdot \log (w_{i}) \\ &\color{steelblue}{\alpha \sim \mathrm{Normal}(178, 100)} \\ &\color{steelblue}{\beta \sim \mathrm{Normal}(0, 100)} \\ &\color{steelblue}{\sigma \sim \mathrm{Exponential}(0.01)} \\ \end{align*} \]
Où \(h_{i}\) est la taille de l’individu \(i\) et \(w_{i}\) le poids de l’individu \(i\). La fonction pour calculer le log en R est simplement log(). Est-ce que vous savez interpréter les résultats ? Indice : faire un plot des données brutes et surimposer les prédictions du modèle…
# on garde seulement les individus ayant moins de 18 ans
d <- open_data(howell) %>% filter(age < 18)
priors <- c(
prior(normal(150, 100), class = Intercept),
prior(normal(0, 10), class = b),
prior(exponential(0.01), class = sigma)
)
mod7 <- brm(
height ~ 1 + weight,
prior = priors,
family = gaussian(),
data = d
) Family: gaussian
Links: mu = identity
Formula: height ~ 1 + weight
Data: d (Number of observations: 192)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept 58.25 1.40 56.04 60.46 1.00 5094 3481
weight 2.72 0.07 2.61 2.83 1.00 5147 3289
Further Distributional Parameters:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma 8.54 0.44 7.87 9.26 1.00 4723 3005
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 5, to = 45, length.out = 1e2) )
# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(
fitted(mod7, newdata = weight.seq, probs = c(0.055, 0.945) )
) %>%
bind_cols(weight.seq)
pred_height <- data.frame(
predict(mod7, newdata = weight.seq, probs = c(0.055, 0.945) )
) %>%
bind_cols(weight.seq)
# on affiche les 6 premières lignes de pred_height
head(pred_height) Estimate Est.Error Q5.5 Q94.5 weight
1 71.84516 8.715536 58.18778 85.59114 5.000000
2 72.99809 8.598611 59.31744 86.57272 5.404040
3 74.13813 8.560623 60.77237 87.66714 5.808081
4 75.33335 8.679844 61.72114 89.19021 6.212121
5 76.21706 8.510558 62.54384 89.86062 6.616162
6 77.26063 8.886644 63.11313 91.46341 7.020202
d %>%
ggplot(aes(x = weight, y = height) ) +
geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
geom_ribbon(
data = pred_height, aes(x = weight, ymin = Q5.5, ymax = Q94.5),
alpha = 0.2, inherit.aes = FALSE
) +
geom_smooth(
data = mu, aes(y = Estimate, ymin = Q5.5, ymax = Q94.5),
stat = "identity", color = "black", alpha = 0.8, size = 1
) Family: gaussian
Links: mu = identity
Formula: height ~ 1 + log(weight)
Data: d (Number of observations: 544)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept -23.59 1.33 -25.72 -21.42 1.00 3790 2568
logweight 47.02 0.38 46.39 47.64 1.00 3809 2666
Further Distributional Parameters:
Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma 5.16 0.15 4.92 5.41 1.00 4752 2789
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 5, to = 65, length.out = 1e2) )
# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(
fitted(mod8, newdata = weight.seq, probs = c(0.055, 0.945) )
) %>%
bind_cols(weight.seq)
pred_height <- data.frame(
predict(mod8, newdata = weight.seq, probs = c(0.055, 0.945) )
) %>%
bind_cols(weight.seq)
# on affiche les 6 premières lignes de pred_height
head(pred_height) Estimate Est.Error Q5.5 Q94.5 weight
1 52.21110 5.260987 43.77279 60.63197 5.000000
2 57.42456 5.190758 49.06822 65.75867 5.606061
3 62.30021 5.158225 54.13302 70.45178 6.212121
4 66.60076 5.284018 58.12409 74.90973 6.818182
5 70.72467 5.160576 62.37382 78.96521 7.424242
6 74.38339 5.218961 66.21511 82.83285 8.030303
d %>%
ggplot(aes(x = weight, y = height) ) +
geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
geom_ribbon(
data = pred_height, aes(x = weight, ymin = Q5.5, ymax = Q94.5),
alpha = 0.2, inherit.aes = FALSE
) +
geom_smooth(
data = mu, aes(y = Estimate, ymin = Q5.5, ymax = Q94.5),
stat = "identity", color = "black", alpha = 0.8, size = 1
)Ladislas Nalborczyk - IMSB2026